clc; clear; close all;
% 在不同 z0 下的分岔图

% 参数设置
mu = 3.5;
k  = 0.16;
x0 = 0.1;
y0 = 0;
z0_range = linspace(-10, 8, 500);

N_transient = 1000;
N_bif = 200;
num_z = length(z0_range);

% 1) 预分配数组大小：每个 z0 要记录 N_bif 个点
totalPoints = num_z * N_bif;
bif_x  = zeros(1, totalPoints);
bif_y  = zeros(1, totalPoints);
bif_z0 = zeros(1, totalPoints);

idx = 1;  % 累计索引
for i = 1:num_z
    z0 = z0_range(i);
    % 瞬态迭代
    state = [x0, y0, z0];
    for t = 1:N_transient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end
    % 记录 N_bif 步
    for t = 1:N_bif
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
        bif_x(idx)  = state(1);
        bif_y(idx)  = state(2);
        bif_z0(idx) = z0;
        idx = idx + 1;
    end
end

% 绘图
figure;
scatter(bif_z0, bif_x, 5, 'b', 'filled'); hold on;
scatter(bif_z0, bif_y, 5, 'r', 'filled');
xlabel('z_0');
ylabel('x, y');
title('变量 x, y 随 z_0 的分岔图');
grid on;
